R Markdown

Lets start by loading in libraries we will need

library("xlsx")
library(ggplot2)
library(forcats)
library(tidyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(priceTools)
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
library(dummies)
## dummies-1.5.6 provided by Decision Patterns

and the data we will be using.

isotope_data <- read.xlsx("edited_community_isotope_data.xlsx",2)

Simple visualization of the data

head(isotope_data)
##   Panel_Site Panel_Number Site_and_Panel              Species_ID     Taxa
## 1         FC            1          FC_01   Alcyonidium polypylum Bryozoan
## 2         FC            1          FC_01   Alcyonidium polypylum Bryozoan
## 3         FC            1          FC_01         Alitta succinea     worm
## 4         FC            1          FC_01   Amphibalanus eburneus Barnacle
## 5         FC            1          FC_01         Bugula neritina Bryozoan
## 6         FC            1          FC_01 Conopeum chesapeakensis Bryozoan
##   pedino_13C_enrichment atom_13C_percent pico_15N_enrichment atom_15N_percent
## 1            -18.113866         1.085858            6.193911        0.3685638
## 2            -20.589898         1.083149            6.303070        0.3686036
## 3            -21.085582         1.082607           12.824082        0.3709834
## 4            -19.677066         1.084148            7.408607        0.3690071
## 5              2.423247         1.108318            6.323871        0.3686112
## 6             -7.060286         1.097948            8.278584        0.3693246
unique(isotope_data$Taxa)
##  [1] "Bryozoan"     "worm"         "Barnacle"     "Bivalve"      "Amphipod"    
##  [6] "Anemone"      "Tunicate"     "brittle star" "sponge"       "crab"        
## [11] "Isopod"

The values contained in the columns pedino_13C_enrichment and pico_15N_enrichment represent FILL IN

The values contained in the columns atom_13C_percent and atom_15N_percent represent FILL IN

Barplots of the Means across each column of interest for Taxa

#barplots by taxa means
par(mfrow=c(2,2))
color_vector_taxa = rev(rainbow(length(unique(isotope_data$Taxa))))
mean_taxa_pedino = tapply(isotope_data$pedino_13C_enrichment, isotope_data$Taxa, mean)
mean_taxa_pedino_percent = tapply(isotope_data$atom_13C_percent,isotope_data$Taxa, mean)
barplot(mean_taxa_pedino,las=2, col = color_vector_taxa, ylab = 'Pedino Enrichment (13C) Mean', main = "Pedino Enrichment Mean vs. Taxa")
barplot(mean_taxa_pedino_percent, las=2, ylim=c(0.0, 1.2), col = color_vector_taxa, ylab = 'Pedino 13C atom% mean', main = "Pedino Atom% vs. Taxa")

mean_taxa_pico = tapply(isotope_data$pico_15N_enrichment, isotope_data$Taxa, mean)
mean_taxa_pico_percent = tapply(isotope_data$atom_15N_percent,isotope_data$Taxa, mean)
barplot(mean_taxa_pico,las=2, col = color_vector_taxa, ylab = "Pico Enrichment (15N) Mean", main = "Pico Enrichment Mean vs. Taxa")
barplot(mean_taxa_pico_percent, las=2, col = color_vector_taxa,ylab = "Pico Atom% (15N)", main = "Pico Atom% Mean vs. Taxa")

Barplots of the Means across variable each column for Site

#barplots by site
par(mfrow=c(2,2))
color_vector_site = rev(terrain.colors(length(unique(isotope_data$Taxa))))
mean_site_pedino = tapply(isotope_data$pedino_13C_enrichment, isotope_data$Panel_Site, mean)
mean_site_pedino_percent = tapply(isotope_data$atom_13C_percent,isotope_data$Panel_Site, mean)
barplot(mean_site_pedino,las=2, col = color_vector_site, ylab = "Pedino Enrichment (13C) Mean", main = "Pedino Enrichment vs Site")
barplot(mean_site_pedino_percent, las=2, ylim=c(0.0, 1.2), col = color_vector_site,ylab = "Pedino Atom% (13C)", main = "Pedino Atom% Mean vs. Site")

mean_site_pico = tapply(isotope_data$pico_15N_enrichment, isotope_data$Panel_Site, mean)
mean_site_pico_percent = tapply(isotope_data$atom_15N_percent,isotope_data$Panel_Site, mean)
barplot(mean_site_pico,las=2, col = color_vector_site, ylab = "Pico Enrichment (15N) Mean", main = "Pico Enrichment Mean vs. Site")
barplot(mean_site_pico_percent, las=2, col = color_vector_site,ylab = "Pico Atom% (15N) Mean", main = "Pico Atom% Mean vs. Site")

Now we will make some boxplots for Taxa and Enrichment

ggplot(data = isotope_data) + 
  geom_boxplot(mapping = aes(x = Taxa, y = pedino_13C_enrichment, color = Taxa)) +  
  labs(x = 'Taxa', y = 'pedino_13C_enrichment', title = 'taxa vs pedino_13C_enrichment') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = isotope_data) + 
  geom_boxplot(mapping = aes(x = Taxa, y = pico_15N_enrichment, color = Taxa)) +  labs(x = 'Taxa', y = 'pico_15N_enrichment', title = 'taxa vs pico_15N_enrichment') + theme(axis.text.x = element_text(angle = 90))

As well as boxplots for Taxa and Atom %

#Pedino
ggplot(data = isotope_data) + 
  geom_boxplot(mapping = aes(x = Taxa, y = atom_13C_percent, color = Taxa)) +  
  labs(x = 'Taxa', y = 'atom 13C %', title = 'taxa vs pedino 13C %') + theme(axis.text.x = element_text(angle = 90))

#pico
ggplot(data = isotope_data) + 
  geom_boxplot(mapping = aes(x = Taxa, y = atom_15N_percent, color = Taxa)) +  
  labs(x = 'Taxa', y = "atom 15N %", title = 'taxa vs pico 15N %') + theme(axis.text.x = element_text(angle = 90))

basic taxa linear regression pedino

pedino_taxa_lm = lm(isotope_data$pedino_13C_enrichment ~ -1 + isotope_data$Taxa)
summary(pedino_taxa_lm)
## 
## Call:
## lm(formula = isotope_data$pedino_13C_enrichment ~ -1 + isotope_data$Taxa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.363  -2.769  -0.437   1.206 182.301 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## isotope_data$TaxaAmphipod        8.132      3.538   2.298   0.0222 *  
## isotope_data$TaxaAnemone       -18.618      7.222  -2.578   0.0104 *  
## isotope_data$TaxaBarnacle      -20.353      1.479 -13.757  < 2e-16 ***
## isotope_data$TaxaBivalve        29.146      5.594   5.210 3.46e-07 ***
## isotope_data$Taxabrittle star   -9.580     12.509  -0.766   0.4444    
## isotope_data$TaxaBryozoan       -3.910      2.553  -1.531   0.1267    
## isotope_data$Taxacrab           -6.921     12.509  -0.553   0.5805    
## isotope_data$TaxaIsopod         -1.754      5.107  -0.343   0.7315    
## isotope_data$Taxasponge         36.403      8.845   4.115 4.96e-05 ***
## isotope_data$TaxaTunicate       -3.953      3.127  -1.264   0.2072    
## isotope_data$Taxaworm           -5.502      2.948  -1.866   0.0630 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.69 on 309 degrees of freedom
## Multiple R-squared:  0.4509, Adjusted R-squared:  0.4313 
## F-statistic: 23.06 on 11 and 309 DF,  p-value: < 2.2e-16

basic taxa linear regression pico

pico_taxa_lm = lm(isotope_data$pico_15N_enrichment ~ -1 + isotope_data$Taxa)
summary(pico_taxa_lm)
## 
## Call:
## lm(formula = isotope_data$pico_15N_enrichment ~ -1 + isotope_data$Taxa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3637 -1.6963 -0.2294  0.8359 24.0512 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## isotope_data$TaxaAmphipod       9.5808     0.5611  17.075  < 2e-16 ***
## isotope_data$TaxaAnemone       10.7529     1.1453   9.388  < 2e-16 ***
## isotope_data$TaxaBarnacle       9.5564     0.2346  40.734  < 2e-16 ***
## isotope_data$TaxaBivalve       10.4804     0.8872  11.813  < 2e-16 ***
## isotope_data$Taxabrittle star   8.3772     1.9838   4.223 3.18e-05 ***
## isotope_data$TaxaBryozoan       8.7510     0.4049  21.611  < 2e-16 ***
## isotope_data$Taxacrab          10.0385     1.9838   5.060 7.20e-07 ***
## isotope_data$TaxaIsopod         5.6950     0.8099   7.032 1.31e-11 ***
## isotope_data$Taxasponge        12.6813     1.4027   9.040  < 2e-16 ***
## isotope_data$TaxaTunicate       9.5369     0.4959  19.230  < 2e-16 ***
## isotope_data$Taxaworm          11.1417     0.4676  23.828  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.805 on 309 degrees of freedom
## Multiple R-squared:  0.924,  Adjusted R-squared:  0.9213 
## F-statistic: 341.6 on 11 and 309 DF,  p-value: < 2.2e-16

Boxplots for enrichment by taxa mean- per panel

FC

FC_isotope <- subset(isotope_data, isotope_data$Panel_Site == "FC")
ggplot(data = FC_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment FC') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = FC_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment FC') + theme(axis.text.x = element_text(angle = 90))

HI

HI_isotope <- subset(isotope_data, isotope_data$Panel_Site == "HI")
ggplot(data = HI_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment HI') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = HI_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment HI') + theme(axis.text.x = element_text(angle = 90))

ID

ID_isotope <- subset(isotope_data, isotope_data$Panel_Site == "ID")
ggplot(data = ID_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment ID') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = ID_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment ID') + theme(axis.text.x = element_text(angle = 90))

IRL_3

IRL3_isotope <- subset(isotope_data, isotope_data$Panel_Site == "IRL3")
ggplot(data = IRL3_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment IRL3') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = IRL3_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment IRL3') + theme(axis.text.x = element_text(angle = 90))

MIM

MIM_isotope <- subset(isotope_data, isotope_data$Panel_Site == "MIM")
ggplot(data = MIM_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment MIM') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = MIM_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment MIM') + theme(axis.text.x = element_text(angle = 90))

MO2

MO2_isotope <- subset(isotope_data, isotope_data$Panel_Site == "MO2")
ggplot(data = MO2_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment MO2') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = MO2_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment MO2') + theme(axis.text.x = element_text(angle = 90))

SCD

SCD_isotope <- subset(isotope_data, isotope_data$Panel_Site == "SCD")
ggplot(data = SCD_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment SCD') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = SCD_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment SCD') + theme(axis.text.x = element_text(angle = 90))

SMS

SMS_isotope <- subset(isotope_data, isotope_data$Panel_Site == "SMS")
ggplot(data = SMS_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment SMS') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = SMS_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment SMS') + theme(axis.text.x = element_text(angle = 90))

WP

WP_isotope <- subset(isotope_data, isotope_data$Panel_Site == "WP")
ggplot(data = WP_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment WP') + theme(axis.text.x = element_text(angle = 90))

ggplot(data = WP_isotope) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment WP') + theme(axis.text.x = element_text(angle = 90))

Now we will create a dummy matrix (isotope_copy) in order to reorder the factor levels of Species ID by Taxa to graph a master plot of the Species and Enrichment values in an organized manner

isotope_copy <- isotope_data
isotope_copy$Species_ID <- factor(isotope_copy$Species_ID)
isotope_copy$Species_ID <- fct_reorder(isotope_copy$Species_ID, isotope_copy$Taxa, .fun = max)

Species versus Pedino Enrichment

ggplot(data = isotope_copy) + 
  geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment master') + theme(axis.text.x = element_text(angle = 90))

Species versus Pico enrichment

ggplot(data = isotope_copy) +
  geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +  
  labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment master') + theme(axis.text.x = element_text(angle = 90))

subsetting the species of algae

pedino_subset <- isotope_data[c('Site_and_Panel',"Species_ID","pedino_13C_enrichment","atom_13C_percent")]
pico_subset <- isotope_data[c('Site_and_Panel',"Species_ID","pico_15N_enrichment", "atom_15N_percent")]

creating community matrixes

comm_pedino <- with(pedino_subset, tapply(pedino_13C_enrichment, list(Site_and_Panel, Species_ID), sum))
comm_pedino <- ifelse(is.na(comm_pedino), 0, comm_pedino)


comm_pico <- with(pico_subset, tapply(pico_15N_enrichment, list(Site_and_Panel, Species_ID), sum))
comm_pico <- ifelse(is.na(comm_pico), 0, comm_pico)

pico pca

pico_pca = rda(comm_pico, scale = TRUE)
pico_pca
## Call: rda(X = comm_pico, scale = TRUE)
## 
##               Inertia Rank
## Total              46     
## Unconstrained      46   25
## Inertia is correlations 
## 
## Eigenvalues for unconstrained axes:
##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
## 7.845 6.618 4.749 3.905 3.617 2.546 2.308 2.109 
## (Showing 8 of 25 unconstrained eigenvalues)
# the eigenvalues sum up to equal the total interia (i.e., total variance in this case)
sum(pico_pca$CA$eig)
## [1] 46
# the ratio of the eigenvalue to the total variance is the amount of 
# variance explained by each PCA axis
round(pico_pca$CA$eig / pico_pca$tot.chi, 2)
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 
## 0.17 0.14 0.10 0.08 0.08 0.06 0.05 0.05 0.04 0.04 0.03 0.03 0.02 0.02 0.02 0.02 
## PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 
## 0.01 0.01 0.01 0.01 0.01 0.00 0.00 0.00 0.00
#We can see from above that the PCA axis 1 captures approximately 17% of the total variance in the community matrix. Let's graph the data to better get a sense of the correlation structure.
plot(pico_pca)

biplot(pico_pca)

ordiplot(pico_pca, display = 'species')
orditorp(pico_pca, display = 'species')

pedino pca

pedino_pca = rda(comm_pedino, scale = TRUE)
pedino_pca
## Call: rda(X = comm_pedino, scale = TRUE)
## 
##               Inertia Rank
## Total              46     
## Unconstrained      46   25
## Inertia is correlations 
## 
## Eigenvalues for unconstrained axes:
##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
## 7.324 6.879 5.468 3.627 3.328 2.813 2.466 2.161 
## (Showing 8 of 25 unconstrained eigenvalues)
par(mfrow= c(1,1))
# the eigenvalues sum up to equal the total interia (i.e., total variance in this case)
sum(pedino_pca$CA$eig)
## [1] 46
# the ratio of the eigenvalue to the total variance is the amount of 
# variance explained by each PCA axis

round(pedino_pca$CA$eig / pedino_pca$tot.chi, 2)
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 
## 0.16 0.15 0.12 0.08 0.07 0.06 0.05 0.05 0.04 0.03 0.03 0.03 0.02 0.02 0.02 0.02 
## PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 
## 0.01 0.01 0.01 0.01 0.00 0.00 0.00 0.00 0.00
#We can see from above that the PCA axis 1 captures approximately 16% of the total variance in the community matrix. Let's graph the data to better get a sense of the correlation structure.
plot(pedino_pca)

biplot(pedino_pca)

ordiplot(pedino_pca, display = 'species')
orditorp(pedino_pca, display = 'species')